% This cript computes tidal power spectra for all the geomagnetic
% stations using two approaches 1) pmtm 2) robustfit

% Get a list of stations

S = dir('/data/backup/mnair/obs_mag_data/update_2011/secular_removed/*.nc');
a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;
t = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30));

for i = 1:length(S),
robstn(i,:) = S(i).name(1:3);
end;
    
[ia,nstn,ib] = intersect(robstn,'CMO','rows')

ncid = netcdf.open(['/data/backup/mnair/obs_mag_data/update_2011/secular_removed/' ... 
    S(nstn).name],'NOWRITE');



X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');
z_data = netcdf.getVar(ncid,Z_ID,0,a);

latitude  = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lat'));
longitude = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lon'));

z_data = double(z_data)/10;
z_data(z_data == 99999.9 ) = NaN;
L = isnan(z_data);
k = find(L==0);

z_data(L) = 0;

netcdf.close(ncid);

if k(end) - k(1) > 500000,
    st = k(1);
    en = k(1) + 500000 - 1;
else
    st = k(1);
    en = k(end);
end;

try
    
% [P, F] = pmtm(z_data(st:en), 4, (en-st) + 1, 1/60);
% 
% LL = 1./(F*3600) > 10 & 1./(F*3600) < 30;
%  
% P = P(LL);
% F = F(LL);
 
display(S(nstn).name);
 
catch me    
end;


%robust fit
%all data
clear k LL

if sum(~L) > 500000,
    L(st+500000:end) = 1;
end;

% %robust fit
% %all data
% 
% spectra_all = robustfit(x(~L,:), z_data(~L));
% 
% % calculate the amplitudes
% 
% spectra_data_all_rob = abs(complex(spectra_all(2:2:end),spectra_all(3:2:end)));


% %least square solution
% 
% spectra_all_ls = x(~L,:)\z_data(~L);
% 
% spectra_data_all_ls = abs(complex(spectra_all_ls(1:2:end),spectra_all_ls(2:2:end)));
% 

% only nightime data

%Get the LT

LL = L | (t-floor(t)  < 18/24 & t-floor(t)  > 6/24)';

% sweeping in frequency for M2 component ROB TOO SLOW
prd = 12.0;


for i = 1:10000,
    
x_m2 = [ cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
    ];  

spectra_night_rob_m2(i,:) = robustfit(x_m2(~LL,:),z_data(~LL));
spectra_data_rob_m2_n(i) = abs(complex(spectra_night_rob_m2(i,2),spectra_night_rob_m2(i,3)));

spectra_night_ls_m2(i,:) = x_m2(~LL,:)\z_data(~LL);
spectra_data_ls_m2_n(i) = abs(complex(spectra_night_ls_m2(i,1),spectra_night_ls_m2(i,2)));

prd = prd + 0.0001;
if mod(i,100) == 0,
    display(i);
end;
end;
periods = 12.415:0.0001:prd;

eval(['save /home/mnair/projects/tides/' S(nstn).name(1:end-3) '_ROB_LS_SWEEP_HD' ...
    ' periods spectra_night_rob_m2 spectra_night_ls_m2 spectra_data_rob_m2_n spectra_data_ls_m2_n st en latitude longitude']);

%clear z_data L LL spectra_day_rob spectra_data_day_rob spectra_data_ls_m2_array spectra_data_night_ls spectra_night_ls spectra_data_all_ls spectra_all_ls spectra_data_rob_m2_n;



